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Abstract 

Cold atoms in optical lattices are described in real space by multi-orbital mean-field Ansatze. In 
this work we consider four typical systems: (i) spinless identical bosons, (ii) spinor identical bosons 
(iii), Bose-Bose mixtures, and (iv) Bose- Fermi mixtures and derive in each case the corresponding 
multi-orbital mean-field energy-functional and working equations. The notions of dressed Wannier 
functions and Wannier spinors are introduced and the equations defining them are presented and 
discussed. The dressed Wannier functions are the set of orthogonal, translationally-equivalent 
orbitals which minimizes the energy of the Hamiltonian including boson-boson (particle-particle) 
interactions. Illustrative examples of dressed Wannier functions are provided for spinless bosonic 
atoms and mixtures in one-dimensional optical lattices. 
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I. INTRODUCTION 



The first experimental realization of the superfTuid to Mott-insulator quantum phase 
transition in cold bosonic atoms trapped by a three-dimensional (3D) optical lattice 
following the theoretical suggestion in |2j], has developed into one of the most active research 
subjects in cold-atom physics nowadays, see the recent reviews {2, 4] and references therein. 
Other experimental studies soon followed, including, e.g., demonstration of the superfluid 
to Mott-insulator transition in effective one- dimensional (ID) optical lattices and obser- 
vation of the so-called Tonks- Girardeau gas thus entering deep into the strongly-interacting 
regime 6] . Very recently, loading of atomic Bose- Fermi mixtures into 3D optical lattices has 
been achieved p], |8| where coherence properties have been measured and studied. 

Trapping of cold atoms in optical lattices (as well as in other trap geometries) occurs in 
real-space. We therefore wish to attack the properties of cold atoms in optical lattices in real 
space. Recently, a multi-orbital best-mean-field ansatz for spinless identical bosonic systems 
has been derived q|. It allowed us to obtain quantitative value of the standard superfluid 
to Mott-insulator transition of weakly-interacting spinless bosons in deep ID optical lattices 
and led us to predict a wealth of quantum phases and excitations of strongly-interacting 
bosons in optical lattices [id ], and novel phenomena associated with fragmentation and 



fermionization of bosons in traps 
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121 ] . Very recently, in a first application of the multi- 



orbital best-mean-field approach to bosonic mixtures, we have described demixing scenarios 
of bosonic mixtures in optical lattices from macroscopic to microscopic length scales 13]. 
The anticipation that the multi-orbital best-mean-field approach would continue to serve as 
a valuable tool for cold-atom systems in optical lattices and other traps serves as the main 
motivation of the present work. To be more specific, the purposes of this paper are: (i) 
to extend the multi-orbital mean-field approach originally developed for spinless identical 
bosons to spinor condensates and to Bose-Fermi atomic mixtures, and to explicitly derive the 
formalism for Bose-Bose atomic mixtures; (ii) to analyze the multi-orbital approach in the 
context of optical lattices, highlighting the role and effects of particle-particle interactions 
and translational symmetry on the properties of the solutions (orbitals). 

The structure of the paper is as follows. In section II we present the 'original' multi-orbital 
mean-field for spinless bosons and discuss some of its properties in the context of bosons in 
optical lattices. Attention is given to the notion of dressed Wannier functions which is intro- 
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duced and demonstrated with the aid of illustrative numerical examples in ID. In section III 
we extend the idea of a multi-orbital mean-field for bosons to spinor condensates. Working 
equations are explicitly derived and discussed in the context of optical lattices. In section 
IV we proliferate the idea of a multi-orbital mean-field for single-species bosons to bosonic 
mixtures. Working equations are explicitly derived. Attention is given to the dressed Wan- 
nier functions which arise due to intra- and inter-species interactions. Illustrative numerical 
examples in ID are provided. In section V we treat Bose-Fermi mixtures, and explicitly 
derive the corresponding working equations of the multi-orbital mean-field approach and 
discuss them in the context of trapped Bose-Fermi mixtures in optical lattices. Finally, in 
section VI we present concluding remarks and some outlook. 

II. MULTI-ORBITAL MEAN-FIELD FOR SPINLESS IDENTICAL BOSONS 
IN OPTICAL LATTICES 

A. Theory 

Our starting point is the many-body Hamiltonian describing N spinless identical bosons 
of a single species in an optical lattice (trap), 



Here, h(r) is the one-body Hamiltonian, containing kinetic and potential energy (the opti- 
cal lattice) terms, is the coordinate of the i-th particle, and AoV^(rj — r^) describes the 
pairwise interaction between the i-th and j-th bosonic atoms where A is a parameter which 
measures the strength of the interparticle interaction. Another assumption, pertinent to 
physically realistic potentials, is that two-body matrix elements associated with Viji — rj) 
(see below) are finite numbers. As is well known from electron-structure theory, even for 
potentials as singular and long-ranged as the Coulomb potential these matrix elements are 
finite quantities. 

As mentioned above, we are going to obtain a real-space mean-field picture of the quantum 
state of cold atoms in the optical lattice. How are we going to achieve that? To this end, 
we attach an orbital to each of the iV atoms. The simplest choice is the Gross-Pitaevskii 



A? 



A' 




(1) 



i=l 



j>i 
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approach, for which all bosons reside in the same orbital, 



tf(r!,r 2 , ...,r N ) 



ri)02(r 2 j 



(2) 



There are, however, many other situations for bosons [9( . Generally, we may take n% bosons 
to reside in one orbital, 0i(r), ri2 bosons to reside in a second orbital, <p2 (r), and so on, 
distributing the N atoms among n or b > 1 orthonormal orbitals. At the other end to the 
Gross-Pitaevskii approach lies the situation where each boson in the optical lattice resides 



in a different orbital, i.e., n or b 
wavef unction 
wavef unction 



N 
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12J. More formally, the multi-orbital mean- field 



x>r N interacting spinless bosons is given by the following single- configuration 

3, 

^(ri, r 2 , . . . , rjv) = S {0i(ri)0 2 (r 2 ) • • • 0jv(rjv)} , (3) 



where S is the symmetrization operator. Note that the Gross-Pitaevskii equation is a specific 
case of Eq. ([3]) when all orbitals are alike, namely n or b = 1. 

To proceed for any set of n or b orbitals and a corresponding set of occupations the 
expectation value of the Hamiltonian ([1]) with the multi-orbital mean-field wavefunction (J3]) 
readily reads [9|, 



E = > //, 



(r)/i(r)0j(r)dr - A, 



rii + 1 



0*(r)0*(r')V r (r - r / )&(r)&(r , )drdr / + 



\ E A n,- J j 0*(r)0*(r')V(r - r') {1 + V rr ,} U*)H*')drdv' 



(4) 



where V rr i permutes the r and r' coordinates of two bosons (particles) appearing to the 
right of it. Note the plus sign preceding the operator V rr > which is due to the (bosonic) 
symmetrization operator. Also, recall here our working assumption that the two-body matrix 
elements are finite quantities. 

The ground state is obtained by minimizing the energy-functional E with respect to its 
arguments which are the number of different orbitals n or b, the set of occupations {n{\ and, 
of course, the shape of the orbitals </>i(r). This results in a set of n or b coupled equations that 



have to be solved self-consistently 



9]. 



Tl rb 



h(r) - A (^ + l)Ji(r) +J2 X o n o l J A r ) + K A r )] \ <M r ) = X)^( r )' 



i 1, . . . , ri or bi 



(5) 



where the direct and exchange potentials are given by 

Jj(r) = j 0*(rOF(r-r%(rO^ 

Kj(v) = J <j>*(r')V(r - r')r rr ^(r')dr' . (6) 

The Lagrange multipliers ^ appearing on the r.h.s. of Eq. (jSJ) are introduced in order 
to ensure orthonormality of the orbitals, j <p*(r)(f)j(r)dr = 5ij, and satisfy the relations 
n i^ij — n-jfiji. It is instructive to mention that, unlike the case of Hartree-Fock equations 
for fermions [lj], the off-diagonal Lagrange multipliers cannot be removed in general; Also 
see subsequent sections and in particular section V dealing with Bose-Fermi mixtures. 

Let us analyze some properties of the energy-functional (j3J) as far as they are needed here, 
assuming V(r — r') to be a repulsive physical interaction. For more details see js), 12]. For 
weakly-interacting atoms, the energy-functional E is minimized when all bosons reside in 
the same orbital (the Gross-Pitaevskii orbital), namely for n or b = 1. When Ao is increased, 
the two-body, mean- field energy terms start to dominate E and the occupations {rii} of 
the bosons would like to become smaller in order to minimize the energy. This leads to 
the occupation of more and more orbitals in the ground state. Eventually, the number of 
orbitals becomes equal to the number of bosons, namely all n.j = 1. In other words, the 
energy functional (j4j) is minimized by a set of n or b = N orthonormal orbitals, a result that 
holds in general trap potentials and for any physical repulsive interaction [q . . 

In what follows we would like to study in some detail the nature and appearance of the 
multi-orbital mean-field solutions in optical lattices. In particular we are interested in how 
the interaction affects, or dresses the interaction-free orbitals. In optical lattices, symmetry 
sets in and it is anticipated that Eq. §5§ possesses solutions which reflect the translational 
and other possible symmetries of the lattice. As usual, periodic boundary conditions are 
assumed, thus we employ a super-cell of a large number N w of potential wells (unit cells). 
To be more specific, we consider an optical lattice with a commensurate filling factor of 
f = jf- bosons per site, and examine solutions to Eq. (jSJ) described by n or b = N w orbitals, 
i.e., Hi = f Vz. In this case and when exploiting the translational symmetry of the lattice, 
we are dealing with a set of N w orthonormal orbitals which may be generated by duplicating 
an orbital </>(r) located in, say, the zeroth unit cell, and translating its replica to all other 
unit cells. Denoting by Rj the set of lattice vectors which generate the entire lattice from 
the zeroth unit cell (set R x = 0), the following relations then hold between the orbitals: 



(r — Rj) = 0j(r) Vz. With these conventions, Eq. ([5]) boils down to, i = 1, . . . , N w : 

N w ~\ N w 

h(r) - (/ + l)AoJ(r) + X °f [ J ( r " R j) + K ( r ~ R i)l \ <K r ) = J2 ^ r + R ~ R ^ 

3 ) 3 

(7) 

What are the solutions of Eq. (JTJ)? If the inter-particle interaction vanishes, i.e., Ao = 0, 
then Eq. (j7j) reduces to, i — 1, . . . , N w : 

N w 

h{r)4>(r) = tH5<K* + R * " R j)> (8) 



which is nothing but the text-book definition of Wannier functions [15|. From this, we 
conclude that the solutions of Eq. (|7j) may be considered as boson-dressed Wannier func- 
tions, namely, Wannier functions that are dressed by the interaction between the bosons. 
From a complementary perspective, the dressed Wannier functions are the set of orthogonal, 
translationally-equivalent orbitals which minimizes the energy of the Hamiltonian including 
boson-boson interaction. 



B. Illustrative examples of boson-dressed Wannier functions 

As mentioned above, the application of the 'original' multi-orbital mean-field to spinless 
bosons in optical lattices has already demonstrated fascinating results, see [h]]. These 
include quantitative determination of the superfluid to Mott-insulator transition in deep ID 
optical lattices and a zoo of Mott-insulator phases and excitations involving higher bands. 
Here, as illustrative numerical examples, we would like to concentrate on the concept of 
boson-dressed Wannier functions, which appear, e.g., in the calculation of self- consistent 



10| ], and study how interaction dresses the bear Wannier functions in 



Mott-insulator phases 
different scenarios. 

We concentrate on ID optical lattices, u(x) = uosm 2 (kx), where k is the wave vector. 
Optical lattice depths, uq, are expressed in terms of the recoil energy, Er = h 2 k 2 /2m, where 
m is the mass of the atoms. The usual contact interaction is employed, XqV(x—x') = Xq5(x- 
x'), where Ao is related to the scattering length and the transverse harmonic confinement 
The strength of the inter-particle interaction is expressed via the dimensionless parameter 
7 = m\o/h 2 n (the ratio of interaction and kinetic energies), where n is the linear density. 
We denote the (commensurate) filling factor by /. 
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Let us begin by dressing the Wannier functions in a deep optical lattice, Uq = 25Er. 
For a commensurate filling of / = 1, N w orbitals, and a weak interaction strength of 7 = 
0.00776002, the dressed Wannier functions are shown in Fig. la. For the above parameters 
as can be seen in the figure, the spatial overlap between two neighboring functions is small. 
For commensurate filling of / = 2, N w orbitals, and a much stronger interaction strength 
7 = 12.7418, the dressed Wannier functions are shown in Fig. lb. Since there are now 
two interacting bosons per orbital, the dressed Wannier functions are wider and lower. In 
addition, two neighboring functions possess a larger spatial overlap, compare Figs, la and 
lb, in spite of the lattice being deep. 

Next, we consider the dressed Wannier functions in a shallow optical lattice, uq = Er. For 
a commensurate filling of / = 1, N w orbitals, and an interaction strength of 7 = 3.491, the 
dressed Wannier functions are shown in Fig. 2. For shallow lattices the (dressed) Wannier 
functions extend beyond the next-nearest neighbor sites and the spatial overlap between two 
neighboring functions becomes substantial, compare Figs. 1 and 2. Consequently, dressing 
due to particle-particle interaction in shallow optical lattices, even for filling factors as low 
as / = 1, is expected to have a larger impact than in deep optical lattices. 

So far, we dressed the Wannier functions of the lowest band in an optical lattice (tech- 
nically, by taking N w orbitals). We would like before concluding this section to present an 
example of two-band bosons-dressed Wannier functions. To this end, we consider an optical 
lattice with a commensurate filling of / = 2, 2N W orbitals, and an interaction strength of 
7 = 12.7609. Minimizing the multi-orbital energy-functional we obtain a set of 2N W equiva- 
lent Wannier functions which are connected one to the other by translation and/or reflection 
symmetry, see Fig. 3. Of course, all 2N W dressed Wannier functions are orthogonal one to 
the others. Since the bosons repel each other quite strongly, the dressed Wannier functions 
(orbitals) would like to reduce their spatial overlap in order to minimize the interaction in 
the multi-orbital energy-functional (j4j). These 2N W dressed Wannier functions are substan- 
tially different from the bear Wannier functions of the first and second band. We remind 
that, in the absence of particle-particle interaction, the bear Wannier functions in a deep ID 
optical lattice resemble the two lowest eigenstates of the Harmonic-oscillator Hamiltonian 
which clearly overlap in each unit cell. 
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III. MULTI- ORBITAL MEAN-FIELD FOR SPINOR IDENTICAL BOSONS 



The above theory of multi-orbital mean-field for spinless identical bosons can be readily 
expanded to treat spinor condensates. This task is carried out below for spin-1 (spinor) 
identical bosons. Our starting point is the following many-body Hamiltonian describing iV 
spin-1 bosons of a single species in a trap (optical lattice), see, e.g., 



N 



N 



H(r 1 ,r 2 ,...,r iV ) = ^U(r l -)l« + J]y(i 



i=i 



j>i 



A l w ®l (j) + A 2 S ^®S^ 



u=x,y,z 



(9) 

Here, h(r) is the one-body Hamiltonian which is assumed to be the same for all three 
components, 1^ is a three-by-three unit matrix operating in hyperfine space on the i- 
th boson. V(Ti — rj) describes the angular- momentum-preserving interparticle interaction 
between the z-th and j-th atoms (for simplicity it is assumed to be a local interaction), where 
Ao and A 2 measure the strength of the interparticle interaction in the spin-independent 
(symmetric) and spin-dependent (asymmetric) parts of (Q. We remark that the contact 
inter-particle interaction V(r$ — r^) = 5(r.; — r^-) is usually taken in the literature 



19. 







20j • For ease of presentation and to avoid displaying multiple tensor-product terms, 1^ 



terms with k ^ i,j are not shown in the Hamiltonian 
usual angular-momentum-1 matrices (ft — 1), 



Finally, the matrices S„ are the 



/n 1 n\ 



1 

1 1 
1 



s _ JL 



-i 

1 -i 
i 



V 1 



(10) 



1 

0-1/ 

In the same manner that we treated spinless bosons, we are going to describe the state 
of the spin-1 condensate in real space. To this end, we again attach a one-body function 
to each of the iV atoms. But this time, for spin-1 bosons, it is a vector of orbitals, or 



a spinor with three components: <fr(r) 



0o(r) 



The normalization condition reads 



f<I>t(r) • <&(r)dr = 1. The simplest choice is the Gross-Pitaevskii approach, for which all 
spin-1 bosons reside in the same spinor, 



*(ri,r 2 ,...,r 



Ni 



*(ri) <g> *(r 2 ) 



Nj 



(11) 



s 



This is the generally employed mean-field ansatz for spin-1 condensates 
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3, Q, 



20]. 



Moreover, in certain cases (e.g., A2 <C Ao) the so-called single-mode approximation has 
been invoked in which the spatial dependence of each component is alike, namely <fr(r) is 



taken as <fr(r) 



|a + | 2 + |a | 2 + \a. 



a - 
a 



l 2 = l. 



, where a±, ao are scalars satisfying the normalization condition 



The most general mean- field is obtained by placing n\ spin-1 bosons in one spinor, $i(r), 
n 2 spin-1 bosons in a second spinor, $2(1"), and so on, distributing the N spin-1 atoms among 
n or b > 1 orthonormal spinors. Note that the individual components of different spinors need 
not be orthogonal to one another. This defines the multi-orbital mean-field ansatz for spinor 
condensates, 



*(r 1; r 2 , . . . , r N ) = S {*i(ri) ® $ 2 (r 2 ) ® • • • ® *jv(rjv)} . 



(12) 



As in the spinless case discussed before, the corresponding Gross-Pitaevskii ansatz (fTTj) is a 
specific case of Eq. ( 1121) when all spinors are alike. 

To proceed, let us derive the multi-spinor energy-functional. For any set of n or b spinors 
and a corresponding set of occupations {rii}, the expectation value of the Hamiltonian (Q 
with the multi-orbital (multi-spinor) mean-field wavefunction (|12p reads 



E 



spin—l 



YlnA J *J(r) • fc(r)*i(r)dr - J j drdr'V{r - r') x 



Ao(*t(')-*i(')) (*I(r')-*i(r')) +A 2 ^(*t( r ). S ^(r)) (<fr|(r') • S^(r')) +^£n i x 
j j dvdv'V{v - v')i\ Q [(*|(r) • *i(r)) (*]( r ') • ^(r')) + (*}(r) ' *J'( r )) (*]( r ' **( r '))] + 
A 2 £ [(*J(r) • S,*,(r)) (*J(r') • S^r')) + (*J(r) • S„*,-(r)) (*J(r') • S^(r'))] }}. (13) 



Once the multi-orbital energy-functional is formulated, the ground state is obtained by 
minimizing E spin ^ 1 with respect to the number of different spinors n or b, the set of occu- 
pations {rii} and, of course, the shape of the spinors $«(r). This results in a set of n or {, 
coupled equations for the spinors <frj(r) [or, 3n or f, coupled equations for their respective 
components +) j(r), 0o,i( r ) an d ^_ i(r)] that have to be solved self-consistently and bear a 
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similar appearance to spinless-boson case, see Eq. (JSJ), 

h(v)l - {n t + l)J,(r) + J> [J,(r) + K,-(r)] ^(r) ^ r). 

3 ) 3 

i = l,...,n orb . (14) 

Here, the direct and exchange matrix-potentials are defined by, 

J,(r) = yV(r - r')dr> A (*j(r') • ^(r')) 1 + A 2 £ (*,V) ■ S^(r')) S„ , (15) 



Kj-(r)*i(r) = J V(r-r')dr' 



A • * 4 (r0) 1 + A 2 X] (*J(rO • S,* 4 (r')) S„ 



The Lagrange multipliers ^ty appearing on the r.h.s. of Eq. (|l4p are introduced in order 
to ensure orthonormality of the spinors, j $*(r) ■ <&j(r)dr = 5y, and satisfy the relations 
— njl^ji- We remind that the individual components of different spinors need not be 
orthogonal to one another. As for the spinless case, the off-diagonal Lagrange multipliers 
cannot be removed in general. 

Analysis of the spin-1 multi-orbital energy-functional E spin _ 1 can be drawn along similar 
lines to that performed in the 'original', spinless multi-orbital case. For instance, we may 
enquire what is the "pathway to fermionization" of a spinor condensate when Ao and A 2 
are enlarged by increasing the transverse confinement in an effective ID trap. Clearly, from 
some interaction strength on the energy-functional (fl3|) is minimized when every spin-1 atom 
resides in a different spinor. Beyond that, detailed investigations are needed to determine 
how the individual components 0±,j(r) and 0cm ( r ) behave in the strongly-interacting limit. 

In the present context of optical lattices and analogously to the spinless case, dressed 
Wannier spinors, *fr(r), are obtained due to particle-particle interaction. By considering 
an optical lattice with a commensurate filling factor of / = spin-1 bosons per site, the 
equation defining the dressed Wannier spinors is obtained from Eq. ffT4j) following the same 
lines employed in section II. The result resembles Eq. (JTj) where spinor quantities substitute 
there spinless ones, i — 1, . . . , N w : 

N w N w 

h(r)l - (/ + l)J(r) + f [J(r " + K(r - R,-)] *(r) = £ A*«*(r + Hi - Rj). 

j ) j 

(16) 

Note that the coupling between the different components 0±(r),0o( r ) arising for A 2 7^ 
results in each component being dressed differently. This finding suggests that distinct 
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dressed Wannier functions may become important for describing the different components 
of spinor condensates in optical lattices. 



IV. MULTI-ORBITAL MEAN-FIELD FOR BOSE-BOSE MIXTURES 

So far, we treated bosonic systems made of one type of atoms. The 'original' multi-orbital 
mean-field in the single-species case has been quite successful for bosons in optical lattices 
and other trap potentials. Very recently, in a first application of the multi-orbital best- 
mean-field approach to bosonic mixtures, we have described demixing scenarios of bosonic 
mixtures in optical lattices from 

the purpose of this section is to provide an explicit and detailed derivation of the Bose-Bose 
multi-orbital formalism (i.e., for a mixture of two kinds of spinless bosons) and to discuss 
some of its general properties in the context of bosonic mixtures in optical lattices. 



A. Theory 

Our starting point is the many-body Hamiltonian describing N = Na + N B spinless 
bosons, Na bosons of type A and N B bosons of type B, in a trap (optical lattice), 

N A N B 

Hab(ti, ...,r NA ,r[,...,r' NB ) = H a + H b + J2J2 X abV ab {Ji ~ r£) ■ (17) 

i p 

Here, H A is the single-species Hamiltonian (pQ) where the corresponding quantities are de- 
noted there by the subscript A, e.g., A^V^r* — r?) denotes the intra-species interaction. 
Similarly, H B is the single-species Hamiltonian of the 5-type bosons where the correspond- 
ing quantities are denoted there by the subscript B. Finally, XabVab^i — r') denotes the 
inter-species interaction between the z-th A boson and the p-th B boson, where Xab measures 
the strength of the inter-species interaction. 

As mentioned above, we are going to obtain a real-space mean-field picture of the quantum 
state of the bosonic mixture in the trap. Similarly to the single-species case, we attach an 
orbital to each of the N = Na + N B atoms. The simplest choice is governed by the so-called 



two-component Gross-Pitaevskii approach, for which all A-t 



Dosons reside in one orbital 



21 



22 



23] 



and all -B-type bosons reside in another orbital, see, e.g., [19l . 

*(n, ...,r NA ,r[,...,r' NB ) = 0(n)0(r 2 ) • • • <f>(r NjL ) x WiWte) ' ' ' Wn b )- (18) 

11 



There can, however, be many other situations for the bosonic mixture, by writing a multi- 
orbital mean- field ansatz for each one of the species in the mixture [3] . Generally, we may 
take Hi ^4-type bosons to reside in one orbital, 0i(r), ri2 A-type bosons to reside in a second 
orbital, </> 2 (r), and so on, distributing the N A atoms among n OT ,f> > 1 orthonormal orbitals. 
Independently, we may take m x .B-type bosons to reside in one orbital, ipi(r), m i -B-type 
bosons to reside in a second orbital, ^(r), and so on, distributing the Nb atoms among 
m or b > 1 orthonormal orbitals. We point out that due to the distinguishability of the A- and 
B-type bosons, no relations are assumed between the 0's and ^'s orbitals which, therefore, 
are a priori allowed to overlap. All together, the multi-orbital mean-field wavefunction 
for the N = Na + Nb bosons in the mixture is given by the multiplication of two single 



configurations 



13|: 



tf(ri, . . . ,r NA ,r[, . . .,r' Ng ) = S {0i(ri)0 2 (r 2 ) ■ • • <p NA (r NA )}xS {^(r'O^O^) • ■ -ipN B (r' NB )} ■ 

(19) 



Note that the two-component Gross-Pitaevskii approach Ji], 21, 22, 23] is a specific case 
of Eq. (Tl9l) when all A-type orbitals are alike and all -B-type orbitals are alike, namely 

Tlorb TTlorb 1 • 

To proceed for any set of n or b and set of m orb orbitals and corresponding sets of occupa- 
tions {rii} and {rrii}, the expectation value of the Hamiltonian (II7|) with the multi-orbital 
mean-field wavefunction ffl9|) readily reads: 

n orb m orb p p 

E AB = E A + E B + J2J2 X ABn imp j j cj)*{r)^*{r')V AB {r - r , )0 ^ (r)^(r / )c^rdr , . (20) 

i p 

Here, E A is the single-species multi-orbital energy-functional, see Eq. (j4j), where the A boson 
quantities are inserted therein, and similarly E B is the single-species energy-functional of 
the 5-type bosons. The last term in the Bose-Bose energy-functional E A b represents the 
interaction between the two species. 

The ground state is obtained by minimizing the energy-functional E A b with respect to 
the numbers of different A- and -B-type orbitals n or b and the sets of occupations {rii} 
and {m p } and, of course, the shape of the orbitals <f>i(v) and i^ p (r). This results in a set of 



12 



n or b + m or b coupled equations that have to be solved self-consistently, 

h A (r) - A A (rii + l)J A)i (r) + ^ X A nj [J A ,j(r) + K Aj (r)} + ^ X AB m q J AB , q (r 




^ orb 

= y]ii A;i j(j)j(r), i = I, . . . ,n orl . 

j 

h B (r) - X B {m p + 1) J BlP (r) + X B m q [J B , q (r) + K B<q (r)} + ^ ^ABnjJ BA j(r 

i j 

= ^2vB, pq ipq(r), p = l,...,m orb . 

q 

Here, J A j(r) and K A j(r) are the intra-species direct and exchange potentials, see Eq. ([6]), 
where the A boson quantities are inserted therein, and similarly J B)Q (r) and K B)Q (r) are 
defined for the 5-type bosons. The inter-species direct-like potentials are given by 

■W r ) = J r q (r')V AB (r-r')^j q (r')dr', 

Jbaj(t) = J (f>*(r')V AB (r - r')^{v')dv' . (22) 

The Lagrange multipliers [i At ij are introduced in order to ensure orthonormality of the A- 
type orbitals, J 4>*(r)(j)j(r)dr = S^, and satisfy the relations Ui\i A ^ = rij\i A ^. Analogously 
and independently \ the Lagrange multipliers /is iP9 are introduced to ensure orthogonality of 
the B orbitals, J ^*(r)^ g (T)dr = 5 pq , and satisfy the relations m v \i Bm = m q \i B ^ qv . As for 
the single-species case, the off-diagonal Lagrange multipliers of both sets cannot be removed 
in general. 

Having at hand the energy-functional E AB we can study the changes in the ground- 
state as the inter-species and/or intra-species interactions are varied. Very recently, as 
mentioned above, in a first application of the Bose-Bose energy-functional (|2"U|) and as a 
specific case-study of the "pathway to fermionization" of bosonic mixtures, we have described 
demixing scenarios of bosonic mixtures in optical lattices 13) . For completeness, we briefly 




discuss below how the physics of Ref. 13] emerges from the Bose-Bose energy-functional 
(f2"U|) . The main advantageous of the multi-orbital mean-field approach is that the orbitals 
{<pi(r)} and {^(r)} are determined self-consistently and thus depend on the intra- and inter- 
species interactions and on the density of each bosonic species, whereas in previous treatment 



of bosonic mixtures in optical lattices, see, e.g., Ref. 24j and references therein, orbitals 



(Wannier functions) whose shapes do not depend on these parameters were employed. 

13 



Consider a given bosonic mixture with some fixed intra-species interaction strengths g A , 
g B and inter-species repulsion gAB- We can minimize the energy-functional (120]) and find 
the ground state, i.e., the set of orbitals {0j(r)}, {^ p (r)} and their respective occupations 
{rii}, {m p } that minimize Eab- For g AB <C g A , Qb the single-species terms E A and E B 
dominate Eab because the inter-species interaction term [the third term on the r.h.s. of 
Eq. (1201) ] is much smaller. Consequently, both species can completely mix and spread all 
over the optical lattice. Moreover, the physics of each species is determined by its own 
parameters almost independently of the other species state. Increasing g AB influences the 
physics of both species. In order to reduce the inter-species interaction energy in ( 1201) . the 
orbitals {0j(r)} and {ip p (r)} reduce their overlap in space. When g AB ^> g A ,g B we expect 
this overlap to vanish, which means that the two species completely separate and occupy 
different regions in space. In between these two extreme cases, there are many intriguing 
possibilities where the self- consistent orbitals of the A- and -B-type bosons 'intermingle' in 
between one another, what results in demixing scenarios on various length scales [13] ]. 

Next, let us discuss the dressing of Wannier functions in a bosonic mixture trapped in an 
optical lattice. To this end, we consider an optical lattice with commensurate filling factors 
of f A and f B , namely n.j = f A \/i and m p = f B Vp. In this case, by exploiting the translational 
symmetry of the lattice we arrive at the final result for the equations defining the dressed 
Wannier functions <p(r) and ^(r) of the bosonic mixture which read, i = 1, . . . , N w and 
p = 1, . . . ,N W : 

Nn, 

h A {r) - \ A (f A + l)J A (r) + Yl [ X ^ A t J ^ r " R ^') + Ka{ * ~ R ^)] + 



N u 



+>^ABfBJAB(r - Rj) j<f){r) = 2_^HA,ij4>{r + R, - R,-), 
h B (r) - X B (f B + 1) J fl (r) + [ X bIb [Jb(t - R,) + K B (r - R q )} + 



N u 



+^ABfAJBA(r - R q ) >ip(r) = ^2/j, BtPq ip(r + Rp - Rq). (23) 

J g 

Obviously, the inter- and intra-species interactions dress together the A- and -B-type-boson 
Wannier functions. In the absence of inter-species interaction, Eq. (I23"j) boils down to two 
independent equations for the single-species dressed Wannier functions of the A bosons and 
of the B bosons, see Eq. (jTJ), whereas in the absence of all interactions it boils down to two 
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independent equations for the text-book Wannier function, Eq. (JHJ). 



131 ]. and 



B. Illustrative examples of Bose-Bose dressed Wannier functions 

As mentioned above, the application of the Bose-Bose multi-orbital energy-functional 
(|20|) to bosonic mixtures in optical lattices has already demonstrated fascinating results in 
connection with demixing, see jl3 |. Here, as illustrative numerical examples, we would like 
to concentrate on the concept of Bose-Bose dressed Wannier functions, which appear, e.g., 
in the self- consistent treatment of mixtures of strongly-interacting Mott-insulators 
demonstrate how inter-species interaction - in addition to the intra-species interactions - 
dresses the bear Wannier functions of each species. 

We consider ID optical lattices and work in dimensionless units in which the one-body 
Hamiltonian of the A- type bosons reads h^x) = — ~-pr + 10 cos 2 (x) where the depth of the 
optical lattice is 20 recoil energies. The same Hamiltonian is taken for the S-type bosons. 
The intra- and inter-species particle-particle interactions, qaYa{x — x'), QbVr(x — x') and 
9abVab(x — x ')i are taken to be the common contact potential, see, e.g., 
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22j,|23j. The 



repulsion strengths, Qa, Qb and gAB, are proportional to the corresponding scattering lengths 



and implicitly include the confining parameters of the transverse directions 16]. /a, /b 
denote the filling factors. In ID one often expresses the strength of the (intra-species) 
interaction in terms of the dimensionless parameter ja which is the ratio of interaction and 
kinetic energies. In the above units one has 'Ja — n 'fj and similarly for the B-type atoms. 

We consider the filling factors = fs = 1, Na = N B = N w orbitals, and the strong 
intra-species interactions Qa = 9b = 10. For the inter-species interaction gAB = 1 the 
dressed Wannier functions of the A- and IB-type bosons are depicted in Fig. 4a. As can be 
seen in the figure the two sets of functions sit one atop the other and cannot be distinguished. 
For the above parameters and stronger inter-species repulsion gAB — 10 the dressed Wannier 
functions of each species move apart, see Fig. 4b; In each site, the two orbitals are pressed 
one against the other and thus narrow, and localize at the borders of the lattice site. Since 
the inter-species repulsion is quite strong, the A- and B-type dressed Wannier functions 
(orbitals) would like to reduce their spatial overlap in order to minimize the interaction 
in the Bose-Bose multi-orbital energy-functional (I2"0l . Finally, it is instructive to contrast 
the dressed Wannier functions of the A- and S-type bosons presented in Fig. 4b and the 
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single-species two-band bosons-dressed Wannier functions depicted in Fig. 3. The latter 
show negative parts in each well, since all 2N W dressed Wannier functions are orthogonal to 
one another. In contrast, there can hardly be seen negative parts of the dressed Wannier 
functions of the Bose-Bose mixture, because dressed Wannier functions belonging to different 
species need not, of course, be orthogonal to one another. 



V. MULTI-ORBITAL MEAN-FIELD FOR BOSE-FERMI MIXTURES 



The final system we wish to study is a mixture made of two kinds of particles of different 
quantum statistics, namely bosons and fermions. Our starting point is the many-body 
Hamiltonian describing a Bose-Fermi mixture of N = Nb + Np interacting particles, Nb 
spinless identical bosons and Np spin-half fermions, in a trap (optical lattice). We assume 
for simplicity that the fermionic one- and two-body parts of the Hamiltonian are spin- 
independent. In this case the Hamiltonian is given by, 

N B Np 

H B f{ti, r NB , r[, r' Np ) = H B + H F + ^ ^W^bf^ - r ' P )- (24) 

i p 

Here, both Hb and Hp are the Hamiltonian of single-species interacting particles, see (OQ), 
where the corresponding quantities are denoted there by the subscripts B and F, respectively. 
^bfVbf^i — r p) denotes the inter-species interaction between the z-th boson and the p-th 
fermion where Xbf measures the strength of the inter-species interaction. 

In the multi-orbital mean-field approach, similarly to the previous sections, we attach 
an orbital to each of the N = Nb + Np atoms. In the simplest case, all bosons reside in 
the same orbital whereas each fermion resides, of course, in a different spin orbital. The 
corresponding mean-field wavefunction, which is symmetric under permutation of any two 
bosons and is anti-symmetric to permutation of any two fermions, is simply given by 

#(n, . . . ,t Nb , r[, . . . , r' Np ) = 0( ri )0(r 2 ) • • • <f){r NB ) x A {M^i)M^) ' ' ' ^(r^)} , (25) 

where A is the anti-symmetrization operator. This standard mean-field appro ach has been 



employed in the literature for Bose-Fermi mixtures in traps, see, e.g., 25l, l26l, l27| . usually 
when the fermionic atoms are taken to be spin-polarized (also see below). It may be termed 
the Gross-Pitaevskii-Hartree-Fock approach for Bose-Fermi mixtures. 
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The most general mean-field for the Bose-Fermi mixture is the one which allows also 
for the bosons to occupy different orbitals. Following the multi-orbital approach in the 
purely-bosonic case, we may take n\ bosons to reside in one orbital, 0i(r), n 2 bosons to 
reside in a second orbital, 02( r ) ; and so on, distributing the Nb bosonic atoms among 
n or b > 1 orthonormal orbitals. All together, the multi-orbital mean-field wavefunction for 
the N = Nb + N F atoms of the Bose-Fermi mixture is the multiplication of two single 
configurations, a general permanent for the bosons and a determinant for the fermions, 

• • -,r NB ,r' 1 , . . .,t' Nf ) = S {0i(r 1 )0 2 (r 2 ) • ■ • (f) NB (r NB )}xA {^(r'J^K) • ■ ■ tp NF (r' NF )} . 

(26) 

Of course, the wavefunction (!26j) possesses the appropriate symmetries with respect to per- 
mutation of any two identical particles. To proceed, we would like to prescribe the energy- 
functional. For the bosons we may take any set of n or b orbitals and a corresponding set of 
occupations {n^}. For the fermions, we remind that the Hamiltonian is spin-independent 
and denote in the following quantities related to spin up and spin down atoms by a and 
(5 superscripts, respectively. The most general single-determinant state is comprised of N F 
spin up fermions and N F = Np — Np spin down fermions, corresponding to an unrestricted 
Hartree-Fock ansatz 14J . It is implicitly assumed that if the fermionic gas is polarized, say 



Np = N and Np = 0, then terms involving spin down fermions drop out in what follows. 
With these conventions, the expectation value of the Hamiltonian with the multi-orbital 
Bose-Fermi mean-field wavefunction (|26|) readily reads: 



Ebf — Es+Ep+y 



J^XBFUi J J 0*(r)(^ a (rO)V BF (r-r')^(r)^ a (rO^r' + {a^/3} 



(27) 



Here, Eb is the single-bosonic-species multi-orbital energy-functional, see Eq. (j3J), where 
the quantities are denoted there with the subscript B, and E F is the standard, unrestricted 
Hartree-Fock energy-functional 



E. 



E 



P L 



Np, 

{^(r)Jh F (r)^(r)dv+^X F J j (^(r)J (^(v')JV F (v - r') {1 - V rr ,} 



x 



<(r)<(r')drdr' 



r/3 



p i 



«(r)JK(r'))V f (r - r')^(r)<(r')drdr'. 

(28) 
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The last term in the Bose-Fermi mixture energy-functional (|2"7|) represents the interaction 
between the bosons and fermions. Finally, it is instructive to remark while examining the 
unrestricted Hartree-Fock energy-functional itself that, for X F V F (r — r') = X F 5(r — r') only 
interactions between spin up and spin down fermions contribute to f[2"Sj) and correspondingly 
to the Bose-Fermi energy-functional f[2"Tl) . 

The ground-state of the Bose-Fermi mixture is obtained by minimizing the energy- 
functional E BF with respect to the number n or b of orbitals in which the bosons reside, 
their occupations {rii}, how many fermions are with spin up and how many are with spin 
down, and, of course, the shape of all bosonic and fermionic orbitals. This results in a set 
of n or b + Np + Np coupled equations that have to be solved self-consistently, 



h B (r) - X B (ni + 1) J B ,i(r) + X B nj [^(r) + K BJ (r)} + X BF 



BF,q\ 



J 



2_j ^BF,q( l 



h F (r) + A, 



WW 



h F (r) + Ai 



g 



£A BF n i J i , Bii (r) J- ^p(r) 
p= l,...,JVf, 

+ £ ^BFnjJ F B,j(r) \ (r) 



V 



(29) 



Here, Jej(r) and K B j(r) are the bosonic direct and exchange potentials, see Eq. (jSJ), where 
the above i? quantities are inserted therein. Analogously, J^(r) and i££ (r) are defined for 
the spin up fermions. The inter-species direct-like potentials are given by, 



^BF,q( r ) 



{^(v')JV BF (v-v')^(r')dr', 
Jfbj(t) = j <P)W)V B F{r-v')^{v')dv'. 
For the spin down quantities, interchange a and (3. 



(30) 
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Let us examine a few properties of the multi-orbital working equations (1291) of the Bose- 
Fermi mixture. As for the previous cases treating systems made of bosons only, the La- 
grange multipliers introduced in order to ensure orthonormality of the bosonic 
orbitals, J 4>*(r)(j)j(r)dr = 5ij, and satisfy the relations riifiij = njfiji. Similarly, the 
off-diagonal Lagrange multipliers cannot be removed in general. The Lagrange multi- 
pliers €p q and e^ q are introduced, respectively, to ensure orthonormality of the spin up 
fermionic orbitals, J (^ipp(r)Jijj q (r)dr = 5 pq , and, independently, of the spin down orbitals, 
J (^(r)) ipP{r)dr = 5 pq . Each of the sets of fermionic Lagrange multipliers can indepen- 
dently be diagonalized, because the corresponding direct and exchange potentials are in- 
variant under unitary transformations of the orbitals. This property is indicated by the 
right-hand arrow in the second and third equations of (T2I?j) . where the Lagrange multipliers 
in their diagonal forms are denoted by e" and e^, respectively. This reminds the situation 
of the standard, unrestricted Hartree-Fock equations |l4|. In solving numerically the multi- 
orbital Bose- Fermi system ( l29|) . it is, of course, preferable to work with the diagonal form 
of the fermionic equations [the last two equations in (1291) ] and thereby reduce the compu- 
tational effort. In this case, the resulting set of fermionic orbitals are generally delocalized 
functions. Having in mind the characterization of (localized) dressed Wannier functions also 
for the fermions, we will keep in the following the Lagrange multipliers matrices in their full, 
non-diagonal forms, e° and e^ q . 

The multi-orbital Bose-Fermi energy-functional (1271) can be employed to study different 
properties of Bose-Fermi mixtures, such as quantum phases, the "pathway to fermionization" 
(of the bosons in the mixture) and, analogously to the case of Bose-Bose mixtures, demixing 
of Bose-Fermi mixtures in optical lattices. Here, we would like to concentrate on optical 
lattices and have a closer look at the dressed Wannier functions. We consider as an example 
an optical lattice with commensurate filling factors of fs bosons and one fermion per site. 
The fermions are taken to be polarized, say, with spin up, that is Np = Np. Starting from 
Eq. and making use of the translational symmetry of the lattice, we arrive at the final 
result for the equations defining the dressed Wannier functions of the bosons and fermions 
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in the mixture, 0(r) and ip a (r), which read, i — 1, . . . , N w and p — 1, . . . , iV w : 

n w 

h B (r) - A b (/b + 1) JbW + £ [A b /b [Jfl(r - R,-) + i^(r - R,)] + Xbf J% f {t - R,)] 

N w 

= ^2fiij4>(r + Ri - Rj), 
Mr) + [A F [J£(r - K q ) - Kf(v - R,)] + A BF / B J FB (r - R,)] I #*(r) = 



>0« 



(31) 



As can be seen, the intra-species and inter-species interactions dress the bosonic Wannier 
functions as well as the fermionic Wannier functions. We remind that for XfVf(t — r') = 
XfS(t — r') the spin-polarized fermions do not interact one with the other; In this case the 
fermionic Wannier functions are dressed only due to the interaction with the bosons. 



VI. CONCLUDING REMARKS 



The multi-orbital mean-field approach has been very successful in describing and predict- 
ing physical phenomena of spinless identical bosons in optical lattices and other traps, and 
very recently, also in a first application to bosonic mixtures in optical lattices. These have 
motivated us in the present work to extend, explicitly derive, and to specifically consider the 
multi-orbital mean-field approach for systems commonly studied in the cold-atom-physics 
literature in general and in the optical-lattice community in particular. Specifically, we have 
described in this paper by multi-orbital mean-field Ansatze (i) spinless identical bosons, (ii) 
spinor identical bosons (iii), Bose-Bose mixtures, and (iv) Bose- Fermi mixtures in real space 
optical lattices. 

In the multi-orbital mean-field approach we attach an orbital to each particle. What are 
the 'rules' for constructing a multi-orbital mean-field ansatz for a given physical system and 
what is the relation between the different orbitals employed? First, orbitals belonging to dif- 
ferent particles bear no a priori relation between them and hence are allowed to overlap (i.e., 
they need not be orthogonal to one another). Second, two identical bosons can sit in either 
the same orbital or in two orthogonal orbitals. Third, each fermion, of course, sits in its own 
spin orbital. Fourth, for spinor identical bosons (and for fermions), it is different two spinors 
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which are orthogonal one to the other; their individual respective components need not be 
orthogonal. Finally, the multi-orbital wavefunction is symmetrized (anti-symmetrized) to 
possess the correct permutational symmetry with respect to interchanging any two identical 
bosons (fermions). 

Having prescribed the multi-orbital wavefunction, the expectation value of the Hamil- 
tonian is calculated, leading to the multi-orbital mean-field energy-functional. What are 
the variational parameters of the energy-functional? For bosons, these are the occupation 
of each orbital (or spinor) and the number of orbitals (spinors). Of course, the orbitals 
(spinors) of all bosons and fermions are themselves variational parameters which are to be 
determined self-consistently. Minimizing the energy-functional with respect to its arguments 
leads to a set of coupled, non-linear, generally integro-differential equations for the orbitals. 
These equations have been explicitly derived for a general inter-particle interaction for the 
above-mentioned systems. 

Particular attention in the context of optical lattices has been given to solutions of the 
multi-orbital equations which possesses translational symmetry. These are the dressed Wan- 
nier functions which appear when interaction sets in. The dressed Wannier functions are 
the set of orthogonal, translationally-equivalent orbitals which minimizes the energy of the 
Hamiltonian including boson-boson (particle-particle) interactions. It has generally been 
shown and specifically demonstrated for spinless bosons and mixtures how intra-species 
interactions, inter-species interaction, and couplings to higher-bands can dramatically alter 
the shape of the Wannier functions in comparison to the bear ones. Moreover, an interesting 
picture is expected for spinor particles, in which different couplings of the individual spinor 
components introduce additional possibilities to dress in a different manner each component 
of the Wannier spinor. 

All the above findings demonstrate the wide potential of the multi-orbital approach for 
cold atoms. We anticipate that the employment of the multi-orbital mean-field Ansatze 
presented here is to produce further valuable understanding and predictions of the physics 
of cold atoms in optical lattices, as well as in other traps. 
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FIG. 1: Boson-dressed Wannier functions in a deep ID optical lattice. Plotted in colors are 
three adjacent dressed Wannier functions <p(x) versus the site index "i". The optical lattice is 
illustrated for guidance by the background sinusoidal curve. Also shown is the spatial density 
p( x ) = TST Si w 4> 2 { x ~ Aj), where X% are the lattice vectors. For convenience, 4>(x) is normalized 
on a segment of length 2ir. The optical lattice is of depth 25Er. Other parameters are: (a) 
Commensurate filling of / = 1, N w = 102 sites and 7 = 0.00776002; (b) Commensurate filling 
of / = 2, N w = 51 sites and 7 = 12.7418. The repulsive interaction between two bosons in a 
site makes the dressed Wannier functions wider and increases the spatial overlap between two 
neighboring functions. See text for more details. 
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FIG. 2: Boson-dressed Wannier functions in a shallow ID optical lattice. Same as in Fig. 1 except 
for the parameters: Optical-lattice depth of Er, commensurate filling of / = 1, N w = 102 sites 
and 7 = 3.491. See text for more details. 
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FIG. 3: Two-band boson-dressed Wannier functions in a deep ID optical lattice. Plotted in colors 
are six adjacent dressed Wannier functions [three pairs of <f>(x) and <j>(—x)] versus the site index "i" . 
The optical lattice is illustrated for guidance by the background sinusoidal curve. Also shown is 
the spatial density p(x) = [(j) 2 (x — Xi) + 4> 2 {—x — JQ)] , where Xi are the lattice vectors. 

For convenience, <fi(x) is normalized on a segment of length 2ir. Parameters are: Optical-lattice 
depth of 25Er, commensurate filling of / = 1, N w = 51 sites and 7 = 12.7609. The two-band 
dressed Wannier functions are substantially different from the bear Wannier functions of the first 
and second band. See text for more details. 
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FIG. 4: Dressed Wannier functions of a strongly-interacting bosonic mixture in a deep ID optical 
lattice. Shown are the dressed Wannier functions (j){x) (scaled by 4) of the A- type bosons in black, 
gray, and blue and the dressed Wannier functions ip{x) (scaled by \) of the £>-type bosons in orange, 
magenta and red. The optical lattice is illustrated for guidance by the background sinusoidal curve. 
The index "i" enumerates lattice maxima. Also shown are the spatial densities of the A bosons 
(in black) p A (x) = ^ ^f™ 4?( x ~ X i) and the B bosons (in red) p B (x) = ^ Y^ w ^ 2 (x - X { ), 
where Xj are the lattice vectors. (j>(x) and ip{x) are normalized along the lattice. The parameters 
are: Commensurate filling factors /a = Jb = 1 5 N w = 16 sites, optical lattice of depth 20Er. The 
inter-species interaction strengths are: (a) g^B = 1. A- and l?-type dressed Wannier functions 
are indistinguishable, (b) gAB = 10. The inter-species repulsion is stronger, causing the A and B 
dressed Wannier functions to reduce their overlap in space and consequently become distinct from 
one another. See text for more details. 
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